function ecog = ecogAr2TransferFunction(ecog, channels, interval, ModOrder,mvartype, value)

% ecogAr2TransferFunction calculates for a desired time series recorded at
% different channels measures of causality. Hence a multivariate
% autoregressive process is fitted to the data leaving a set of ar
% parameters (N parameters: NChannels x NChannels x ModOrder) which are
% sufficient to reconstruct the whole dataset. This set of parameters is
% transformed into the frequency domain by multiplying each ar parameter of
% a given ModOrder by exp(-i*2*pi*hz*j*samplinginterval). Per default the
% frequencies range between 1 and 200 Hz. 
% Internally the ecogGrangerTest function is called, which shifts the data
% set channels independently and randomly modulo length of each channel. #
% Function assumes that the data are epoched between -1 and 2 sec around
% the stimulus trigger and trials are concatenated in the third dimension.
% 
% 
% INPUT:    (1) ecog structure containing the variables data and srate
%           (2) channels: subset of recorded channels            
%           (3) interval: Start and End point of interval in which
%                         causality is assumed
%           (4) ModOrder: determines the Order of the ar model
%           (5) mvartype: bayesian ('arfit') or yule-walker ('yule') approach
%           (6) value: this input refers to the theta or the standardized
%                      gamma value. Note that these values have nothing to do
%                      with a special frequency band.
%
%
% OUTPUT:   (1) ecog.H: Transfer function
%           (2) ecog.theta2: abs(ecog.H)^2
%           (3) ecog.gamma2: standardized ecog.theta2 parameter by all
%                            remaining ar parameters of a given channel
%
%
% OUTPUT of internally called ecogGrangerTest:
%           (1) ecog.thetaCI: upper ci cutting 5 % of the distribution
%           (2) ecog.gammaCI: upper ci cutting 5 % of the distribution
%
%
%
% EXAMPLE:  
%           ecog.srate     = 1000;
%           t = (0:ecog.srate-1) * 1/ecog.srate;
%           ecog.data(1,:) = 0.7*sin(2*pi*50*t);
%           ecog.data(2,:) = circshift(ecog.data(1,:)',10);
%           ecog.data(3,:) = circshift(ecog.data(2,:)',10);
%           ecog = ecogAr2TransferFunction(ecog,[1:3], [0 1], 10,'arfit','theta');
%
%
%
% SD wrote it: 2010/08/05


addpath('/data2/duerschm/tmpskripts/arfit')
if isfield(ecog,'H') == 1;
    ecog = rmfield(ecog,'H');
    ecog = rmfield(ecog,'gamma2');
    ecog = rmfield(ecog,'gammaTest');
    ecog = rmfield(ecog,'gammaCi');
end


clc;
% =========================================================================
% general stuff
% =========================================================================
x           = mean(ecog.data,3);
% x = ecog.data(:,:,1);
x           = zscore(x,[],2);
stimOnset   = size(x,2)/3;
StartIdx    = round(stimOnset + stimOnset * interval(1));
EndIdx      = round(stimOnset + stimOnset * interval(2));
dtset       = x(channels, StartIdx:EndIdx);
center_freq = 1:200;


% =========================================================================
% choice between the bayesian and yule-walker approach for fitting a
% multichannel ar model
% =========================================================================
switch mvartype
    case 'arfit'
        [w,A] = arfit(dtset',ModOrder,ModOrder);
    case 'Yule'
        addpath('/data2/duerschm/tmpskripts/');
        A           = mvar(dtset',ModOrder,1);
end
  
% =========================================================================
% concatenating the ar coeffs for each model order in the third dimension
% =========================================================================
A = reshape(A,size(dtset,1),size(dtset,1),size(A,2)/size(A,1));

% =========================================================================
% Transformation of the ar model in the frequency domain
% =========================================================================
freq        = 1;
for hz = center_freq;
    % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    % determining the Transferfunction based on A
    % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    for j            = 1:size(A,3);
        X(:,:,j)     = ((A(:,:,j) * exp(-1i*(2*pi*hz*j*(1/ecog.srate)))));
    end
    Y = inv((sum(X,3) - eye(size(X,1),size(X,1))));
    ecog.H(:,:,freq) = Y;
    freq             = freq + 1;
end

% =========================================================================
% calculation theta_ij^2 as the squared abs value of H_ij and gamma_ij^2
% =========================================================================
x = 1;
ecog.theta2 = zeros(length(channels)^2,length(center_freq));
ecog.gamma2 = zeros(size(ecog.theta2));
for l = 1:size(ecog.H,1)
for m = 1:size(ecog.H,2)
    ecog.theta2(x,:) = squeeze(abs(ecog.H(m,l,:)).^2)';
    ecog.gamma2(x,:) = (ecog.theta2(x,:)./sum(squeeze(abs(ecog.H(:,l,:)).^2),1))';
    x = x + 1;
end
end



% =========================================================================
% calculating null distribution.
% =========================================================================
ecog = ecogGrangerTest(ecog, dtset, 1000,ModOrder,mvartype);


% =========================================================================
% pure plotting
% =========================================================================
channel1 = repmat(channels', 1, length(channels));
channel2 = repmat(channels, length(channels), 1);
switch value
    case 'theta'
        
        figure;
        for k = 1:size(ecog.theta2,1);
            subplot(length(channels), length(channels), k);
            bar(ecog.theta2(k,:));
            hold on;
            plot(ecog.thetaCi(k,:),'r','Linewidth',3);
            set(gca,'ylim',[0 max(max(ecog.theta2))],'xlim',[0 center_freq(end)]);
            title(['ch ',num2str(channel2(k)),' ==> ch ',num2str(channel1(k)),''])
        end
        
        
    case 'gamma'
        
        figure;
        for k = 1:size(ecog.gamma2,1);
            subplot(length(channels), length(channels), k);
            bar(ecog.gamma2(k,:)); 
            hold on;
            plot(ecog.gammaCi(k,:),'r','Linewidth',3);
            set(gca,'ylim',[0 max(max(ecog.gammaCi))],'xlim',[0 center_freq(end)]);
            title(['ch ',num2str(channel2(k)),' ==> ch ',num2str(channel1(k)),''])
        end
end
        



% =========================================================================
% END OF FUNCTION
% =========================================================================